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We use high-temperature series expansions to obtain thermodynamic properties of the quantum 
compass model, and to investigate the phase transition on the square and simple cubic lattices. On 
the square lattice we obtain evidence for a phase transition, consistent with recent Monte Carlo 
results. On the simple cubic lattice the same procedure provides no sign of a transition, and we 
conjecture that there is no finite temperature transition in this case. 

PACS numbers: PACS Indices: 05.30.-d,75.10.-b,75.10.Jm,75.30.Cr,75.30.Kz 
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I. INTRODUCTION 



Quantum compass models are spin models in which 
the nearest-neighbour exchange coupling has the form 
JaSfS'j where a(= x,y,z) depends on the direction of 
the particular link or bond. This then implies a coupling 
between the spin space and the physical space of the lat- 
tice. Such models were first introduced, and have been 
regularly employed, to describe orbital ordering in var- 
ious transition metal compounds ([U-IH, an0 - references 
therein). 

Such models also have applicability to models of p + 
ip superconducting arrays [a, [6( and it has been argued 
[7j that such arrays can provide fault-tolerant qubits for 
quantum information systems. 

Compass models can be defined in various ways, de- 
pending on the underlying lattice. Exact solutions have 
been obtained for a 1-dimensional alternating (xx),(zz) 
model || and for a 2- leg ladder Q. A remarkable solu- 
tion has also been found for the honeycomb lattice with 
(xx) , (yy) and (zz) couplings along the three independent 
lattice directions [l(|. As far as we are aware, no other 
exact solutions exist. 

In the present paper we consider the spin-1/2 quantum 
compass model on the simple cubic lattice, with Hamil- 
tonian 



H = J;, 



(x) 

£• 

<ij> 



(If) 
<ik> 



<il> 



(1) 



where the erf are Pauli operators, and the sums are, re- 
spectively, over lattice bonds along the x,y,z directions. 
We will also consider the square lattice version, where 
the last term in {I} is omitted. 

As is well known [y] , this model possesses a number of 
unusual gauge- like symmetries. As a consequence each 
energy state has a macroscopic degeneracy and, con- 
sequently, there is no conventionally ordered magnetic 
phase at any temperature. However it has been pointed 
out that a state of orientational or 'nematic' order is pos- 
sible, in which the nearest neighbour bonds of lowest en- 
ergy lie predominantly along a specific lattice direction. 
In the isotropic case of equal interactions ( J x = J y = J z ) 



this represents a spontaneous symmetry breaking. Con- 
sequently there may be a critical point at a temperature 
T c , above which the system is disordered, with no pre- 
ferred direction. 

Recent quantum Monte Carlo studies of the isotropic 
2-dimcnsional model have found strong evidence 

for a finite temperature critical point with kT c / J = 0.234 
(in our units) with a critical exponent v ~ 0.97, consis- 



tent with 2D Ising behaviour. The same authors [13 1 
also identified a transition in the corresponding classical 
model, but we do not consider the classical case in the 
present work. As far as we are aware, no investigation of 
the occurrence of such a critical point, or its value, has 
been reported in the 3-dimensional case. 

The goal of the present work is to attempt to an- 
swer this question. We employ the method of high- 
temperature series expansions, which has proven success- 
ful in the past [l4| in obtaining accurate values for critical 
temperatures and exponents in a wide variety of classi- 
cal and quantum models. The basic idea is to expand 
the Boltzmann factor e~@ H in the partition function in 
powers of f3 — 1 / kT 



Z = Trie' 011 } 



Z. — J J*] 



(2) 



r=0 



The coefficients in this series can be evaluated in a num- 
ber of (related) ways. We use a linked cluster approach 
[TEj in which InZ is evaluated, as a series in j3, on a se- 
quence of finite connected clusters of increasing size, and 
the cluster contributions are combined appropriately to 
give the bulk free energy in the form 



BF = —InZ 

N 



= ln2 + a>r(Jxi Jyt Jz)P r 



(3) 



with the a r being multinomial expressions of degree r in 
the J's. From this one can immediately obtain a corre- 
sponding series for the specific heat. 
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However the specific heat has, in most cases, only a 
weak singularity and is not well suited to estimation of 
critical properties. Including an external field which cou- 
ples to the order parameter D, 



H = He] 



hD 



(4) 



where we now write the original Hamiltonian (JTJ as Ho, 
allows calculation of a high temperature series for a gen- 
eralized 'susceptibility' 



(5) 



The order parameter D was introduced [13J for the 2- 
dimensional model as 

0) (y) 
D 2d = J*J2 ~ J vJ2 *M ( 6 ) 

<ij> <ik> 

i.e. the difference between the energy of the x and y 
bonds. We generalize this for the 3-dimensional model 
to 

(2) (*) (y) 

D 3d = 2 J z £ of of - J x g of - J a g ctfo* (7) 

<ij> <ik> 

Normally the calculation of the susceptibility would 
be somewhat involved, since Ho and D do not commute. 
However, in the present model , we can simply combine 
the two terms into a Hamiltonian of the original form ((TJ) , 
with J x -> J x (l - h), J y ->■ J y (l - h), J z -> J z (l + 2h) 
and use the expression in ([3]) to obtain 



oo 



r=2 



(8) 



where the c r are again multinomials of degree r in the 
J's. The susceptibility series is expected to show a strong 
divergence at the critical point and hence should be more 
amenable to analysis. 

Another quantity which is expected to show a strong 
divergence is the fluctuation in the order parameter 



Q =< D 2 > - < D > 2 



(9) 



For the classical model this quantity is identical to \i DU t 
this is not the case for the quantum model. 

In the following sections we will present the series and 
our analysis for the 2-d case (Section [TTJ) and 3-d case 
fSection lllip . Our conclusions are summarized in Section 

m 



II. THE SQUARE LATTICE 

To test the effectiveness of the high-temperature series 
approach for the present model, we first investi gat e the 
square lattice case, where previous results exist [lll - [l3j . 

We use a linked cluster method [15| based on connected 
clusters ('graphs'), To obtain a series for (lnZ)/N correct 
to order /? , as we have done, requires the enumeration 
of clusters with up to 12 bonds. It is a special feature of 
this model that each bond must be used an even number 
of times to give a nonzero trace. There are 4423 topologi- 
cally distinct clusters with 12 or fewer bonds, embeddable 
on the square lattice. This gives rise to 751663 distinct 
graphs with 2 bond types (x and y). However the vast 
majority of these do not contribute, and the final irre- 
ducible list of contributing graphs numbers 60127. We 
give below the leading terms in the partition function 
series 



UnZ = ln2 + l(x 2 +y 2 )f3 2 ~^(x i + 8x 2 y 2 +y i )P i 

+ ^((x 6 + y 6 ) + 30(:rV + ^V))/? 6 - ^(17(^ 8 + V 8 ) + 1376(x 6 y 2 + x 2 y 6 ) + 4344x 4 y 4 )/3 8 
45 2521) 

' ' (31(a; 10 + y 10 ) + 5570(x 8 y 2 +a; 2 2/ 8 )+40500(a;V + a; 4 y 6 ))/3 10 



14175 



' (691(x 12 + y 12 ) + 241800(a; 1 V + x 2 y w ) + 3426402(a; 8 y 4 + x i y 8 ) + 7679480x 6 y 6 )(3 12 ■ ■ ■ (10) 



935550 



where x = J x .y = J y . Note that only even powers of /3 occur. This is a feature of all series for this model. 
From this result we can obtain the susceptibility 

XIP = (x 2 + V 2 ) + y 4 ) - &x 2 y 2 )? 2 + ^((x 6 + rf) 2(xV + zV))/3 4 

--i-(119(a; 8 + y 8 ) + 1376(xV + x 2 y 6 ) - 4344xV)/3 6 

O-l-O 

' (2790(x 10 + y w ) + 144820(x 8 ?/ 2 + a; 2 y 8 ) - 243000(a; 6 y 4 + x 4 y 6 ))/3 8 • • • (11) 



14175 
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TABLE I: Series coefficients for the isotropic 2d Compass Model 



p 


1 i r-7 

jflnZ 


x/WJ 2 ) 


Q 





0.693147180560D+00 


0.200000000000D+01 


0.200000000000D+01 


2 


0.100000000000D+01 


0.666666666666D+00 


0.600000000000D+01 


4 


-0.833333333333D+00 


-0.133333333333D+01 


-0.120000000000D+02 


6 


0.137777777778D+01 


0.429841269841D+01 


0.307111111111D+02 


8 


-0.282936507937D+01 


-0.132382287013D+03 


-0.849650793664D+02 


10 


0.650455026455D+01 


0.421280743947D+02 


0.245585890085D+03 


12 


-0.160518048207D+02 


-0.132382287013D+03 


-0.730131587978D+03 


14 


0.416028785294D+02 


0.418143457749D+03 


0.221416777392D+04 


16 


-0.111781974764D+03 


-0.132765859211D+04 


-0.681447081860D+04 


18 


0.308758184039D+03 


0.423612006027D+04 


0.212139054331D+05 


20 


-0.871688240896D+03 


-0.135761546397D+05 


-0.666457337381D+05 


22 


0.250500394206D+04 


0.436833196403D+05 


0.210941032291D+06 


24 


-0.730521400959D+04 







The higher order terms were evaluated numerically. In 
Table U we show the full series for the isotropic case J x — 
J y . The expansion variable is K — ft J . 

We have attempted to analyse these series using stan- 
dard Pade approximant methods. Our discussion is con- 
fined to the x series, as this (together, possibly, with Q) 
is expected to have a strong singular behaviour at the 
critical point. The first point to make about the series 
in /3 2 is the regular alternation in sign. This reflects the 
presence of a dominant singularity on the negative 1 
axis (i.e. the imaginary /3 axis). In fact there appears 
to be a whole string of such imaginary poles in the Dlog 
Pade approximants. This, in itself, is not so unusual. 
Recall that the exact result for the ID Ising model has 
poles at j3J = ±i(n + l/2)n. 

However these interfering singularities mask the ex- 
pected physical singularity on the real positive axis. 
One possible strategy to overcome this is to use an Euler 
transformation of the form y — x/(l + ax), (x = K 2 ), 
which has the effect of compressing the positive real axis 
and expanding the region — X/a < x < of the negative 
real axis. The use of such transformations is well known 
in the field of critical phenomena, as are the possible pit- 
falls. 

To provide the reader with some insight into the ana- 
lytic structure of the x series we discuss the location of 
poles of Dlog Pade approximants to the series for /3J 2 /x 
before and after the Euler transformation (with a=2.0). 
The original series in x = /3 2 has very consistent poles 
at x ~ —0.28, —0.32, —0.46, with less consistent poles 
much further from the origin. The transformed series 
shows images of these at y = —0.65, —0.9 as well as 
poles on the positive real axis at y = 0.47, 0.54. The last 
of these corresponds to a large negative value x ~ —6.8, 
whereas y = 0.47 corresponds to x — 7.8, or a physical 



TABLE II: Poles and residues (in brackets) in the variable K 2 
for [N/D] Dlog Pade approximants to the quantity x/(PJ 2 ) 
for the 2D compass model, after an Euler transform with a = 
2.0. (Asterisks denote a complex pair of poles in the physical 
region) 



D\N 


3 4 5 6 7 


3 


0.439(0.191) * 0.457(0.288) * 


4 


0.470(0.460) 0.474(0.533) 0.472(0.509) 0.475(0.567) 


5 


0.473(0.529) 0.473(0.514) 0.473(0.575) 


6 


0.473(0.512) 0.473(0.524) 


7 


0.474(0.546) 



critical value kT c / J ~ 0.34. In Table ILTl we show the es- 
timates of y c and the exponent 7 at various orders. As 
can be seen, these are quite consistent at y c ~ 0.473 and 
7 ~ 0.52. However, this critical temperature is much 
higher than the Monte Carlo estimate 0.234 and the cor- 
responding exponent is much lower than the expected 
Ising value of 1.75. Therefore we can only conclude that, 
while the Dlog Pade analysis provides evidence for a 
physical critical point, the numerical estimates cannot 
be taken with any confidence. We comment further on 
this in the conclusions. 

An alternative approach to analysing our series data is 
to evaluate the susceptibility itself at temperatures above 
T c , using Pade approximants, and to plot the inverse 
susceptibility x~ l versus T. In Figure Q] we plot both 
j3/x, obtained directly from the series (fTTj) . and 1/x ver- 
sus temperature. Both curves clearly approach zero at 
T c ~ 0.25, a value consistent with the Monte Carlo es- 
timates [HI, EH , and considerably below our Dlog Pade 
results. It is not possible to obtain accurste exponent es- 
timates from this procedure, but if we fit our data points 
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with a simple form 1/x = a(T — T c ) 7 together with the 
Monte Carlo critical point T c — 0.234 we obtain 7 ~ 1.3, 
which is at least a good deal closer to the expected Ising 
value. 

Thus we conclude that the series approach does con- 
firm the existence of a finite temperature critical point in 
the isotropic 2D model, and corroborates the presumably 
more accurate Monte Carlo results. 



III. THE SIMPLE CUBIC LATTICE 



FIG. 1: Estimated values of the inverse susceptibility /3/x 
(circles) and x _1 (squares) as functions of temperature T for 
the 2D quantum compass model (J=l). 



We now turn to the 3-dimensional model, where no 
previous results exist. We use the same approach as for 
the 2D case, and compute series for the same quantities. 
The leading terms of the series for InZ are 



UnZ = ln2 +^(x 2 + y 2 = z 2 )/3 2 - ^(x 4 + y 4 + z 4 + 8(x 2 y 2 + x 2 z 2 + y 2 z 2 ))f3 4 

+ l.( x * + y 6 + z 6 + m{x A y 2 + x 2 y 4 + X A Z 2 + ^4 + y i J + y 2 ^ + ^^0* 

45 

--i-(17(a; 8 + y s + x s ) + 1376(ssV + x 6 z 2 + y 6 z 2 + y 6 x 2 + z 6 x 2 + z 6 y 2 ) + 4344(»V + x 4 z 4 + y 4 z 4 ) 
2520 

+ 14176(xV^ 2 + y 4 x 2 z 2 + z 4 x 2 y 2 ))p % 

+ ^^(31(x 10 + y w + z w ) + 5570(aV + x s z 2 + y s x 2 + y s z 2 + z s x 2 + z 8 y 2 )) 
14175 

+40500(a;V + y 6 x 4 + x 6 z 4 + z 6 x 4 + y 6 z 4 + z 6 y 4 ) + 120320(x 6 y 2 z 2 + y 6 x 2 z 2 + z 6 x 2 y 2 ) 
+297200(x 4 y 4 z 2 + x 4 y 2 z 4 + x 2 y 4 z 4 ))p w + ■■■ (12) 



where x = J x ,y = J y , z = J z . 

The susceptibility corresponding to the order parame- 
ter D 3d (equation 10) can be obtained by the substitu- 
tion J x -» J x (l — A), J y -> J y (l - A), J z J z (l + 2A) 
in (HJ). This definition, of course, introduces a preferred 



direction z. However in the isotropic limit the resulting 
series is unaffected by this. 

We have evaluated the series numerically, up to order 
/3 20 , and the coefficients are shown in Table Hill 



As for the 2D case, the series are dominated by sin- 
gularities on the negative /3 2 axis. However, in contrast 



to the 2D case, Euler transformations yield no indica- 
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TABLE III: Series coefficients for the isotropic 3D Compass Model 



p 






Q 





0.693147180560D+00 


0.600000000000D+01 


0.600000000000D+01 


2 


0.150000000000D+01 


- 0.600000000000D+01 


0.180000000000D+02 


4 


-0.225000000000D+01 


0.200000000000D+02 


-0.760000000000D+02 


6 


0.673333333333D+01 


-0.810476190476D+02 


0.377466666667D+03 


8 


-0.253440476190D+02 


0.367149206349D+03 


-0.200491428396D+04 


10 


0.107871111111D+03 


-0.178751576719D+04 


0.110779369318D+05 


12 


-0.496475097002D+03 


0.915575874989D+04 


-0.628679773726D+05 


14 


0.241283362972D+04 


-0.486884786086D+05 


0.363814737295D+06 


16 


-0.122062709687D+05 


0.266451000791D+06 


-0.213714381541D+07 


18 


0.636830117143D+05 




0.127041779973D+08 


20 


-0.340463327677D+06 







0.3 



0.2 
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FIG. 2: Estimated values of the inverse susceptibility /3/\ 
(circles) and x _1 (squares) as functions of temperature T for 
the 3D quantum compass model (J=l). 



tion of any singularity for real positive /3 2 , and thus no 
indication of a physical critical point. 

To test this further we have employed the same strat- 
egy as in the previous section, by evaluating x itself at 
high temperatures, where Pade approximants to the se- 
ries are well converged, and plotting versus T. The 
results are shown in Figure [2] 

We note that the f3/x points are monotonically increas- 
ing, unlike the results for the 2D case - Figure [T] This 
indicates that \ is increasing less rapidly than 1/T. The 
1/x values do not indicate a transition at any finite T 



but, within the numerical uncertainties, are consistent 
with a transition at T = 0. 



IV. DISCUSSION 

The question of the existence of a thermodynamic 
phase transition in the quantum compass model on vari- 
ous lattices is of fundamental importance. 

The present work is, to our knowledge, the first at- 
tempt to address this problem using the technique of 
high-temperature series expansions, a standard method 
in other contexts. 

The series indicate that the analytic structure of ther- 
modynamic functions for these models is dominated by 
singularities on the imaginary j3 axis (0 = 1/kT). This 
is perhaps a reflection of the peculiar '1-dimensional' na- 
ture of the couplings in the model. 

Our results for the square lattice are consistent with, 
albeit less precise than, recent Monte Carlo results [l3| . 
This demonstrates that the high-T series method does 
in fact work. However for the cubic lattice we find no 
signature of a critical point at finite T, and conjecture 
that there is no such critical point. At first glance this 
appears surprising, since the normal expectation is that 
the ordered phase will be more robust, and hence T c will 
increase, with increasing dimension. In the case of a sim- 
ple antifcrromagnet, for instance, the bond interactions 
in different directions can be satisfied simultaneously, and 
reinforce each other, so that the tendency to order in- 
creases with higher dimension. In the present case, how- 
ever, the bond interactions in different directions pull 
different ways, and compete with each other, so that the 
tendency to order decreases with higher dimensions. In 
one dimension, the 'nematic' order parameter is non-zero 
at all finite temperatures; in two dimensions Did is only 
non-zero at low temperatures; and in three dimensions 
it appears that D^d is actually zero at all finite temper- 
atures. It has also been pointed out 0, that in this 
model thermal fluctuations in fact become larger with 
increasing dimension. 



6 



The series have proved difficult to analyze, because of 
the complex singularities, and gave rather poor estimates 
of the critical parameters in two dimensions. A closer in- 
vestigation of the nature of these singularities may lead 
to more precise estimates of the critical parameters; or 
else higher-order series coefficients might be necessary. It 
is worth noting that the model has also proved difficult to 
analyze using finite- size scaling and Monte Carlo meth- 
ods. An early Monte Carlo calculation [TTj] on lattices 
of up to 20 x 20 sites with periodic boundary conditions 
also gave a critical point about 36% too high. Wenzel 



et al. [13] showed that the use of special 'screw peri- 
odic' boundary conditions on lattices up to 42 x 42 was 
required to produce the estimates quoted earlier. 



Acknowledgments 

We are grateful for the computing resources provided 
by the Australian Partnership for Advanced Computing 
(APAC) National Facility. 



[1] K.I. Kugel and D.I. Khomskii, Sov. Phys. Usp. 25, 231 
(1982). 

[2] D.I. Khomskii and M.V. Mostovoy, J. Phys. A36, 9197 

(2003) : M.V. Mostovoy and D.I. Khomskii, Phys. Rev. 
Lett. 92, 167201 (2004). 

[3] J. van der Brink, New J. Phys. 6, 201 (2004). 
[4] G. Jackeli and G. Khaliullin, Phys. Rev. Lett. 102, 
017205 (2009). 

[5] C. Xu and J.E. Moore, Phys. Rev. Lett. 93, 047003 

(2004) . 

[6] Z. Nussinov and E. Fradkin, Phys. Rev. B71, 195120 

(2005) . 

[7] B. Doucot, M.V. Feigelman, L.B. Ioffe and A.S. Iosele- 

vich, Phys. Rev. B71, 024505 (2005). 
[8] W. Brzezicki, J. Dziarmaga and A.M. Oles, Phys. Rev. 

B75, 134415 (2007). 



[9] W. Brzezicki and A.M. Oles, Phys. Rev. B80, 014405 
(2009). 

[10] A. Kitaev, Ann. Phys. 321, 2 (2006). 
[11] T. Tanaka and S. Ishihara, Phys. Rev. Letts. 98, 256402 
(2007). 

[12] S. Wenzel and W. Janke, Phys. Rev. B78, 064402 (2008). 
[13] S. Wenzel, W. Janke and A. Lauchli, Phys. Rev. E81, 
066702 (2010) 

[14] C. Domb and M.S. Green (eds.), Phase Transitions and 
Critical Phenomena, Vol. 3 (Academic, New York, 1974). 

[15] J. Oitmaa, C. Hamer and W. Zheng, Series Expansion 
Methods for Strongly Interacting Lattice Models (Cam- 
bridge University Press, 2006). 

[16] A. Mishra, M. Ma, F-C. Zhang, S. Guertler, L-H. Tang 
and S. Wan, Phys. Rev. Letts. 93, 207201 (2004). 



